% Script to predict local magnetic declination variation from 
% local F, and remote XYZ at two obs

% PBC to TAL is 192 km 
% TAL to GHC is 140 km 
% PBC to GHC is 252 km

%% load data.

basedir = '/Users/manojnair/data/obs_mag_data/Canada/Absolute_netCDF/';
start_index = 0; % start from the beginning
ndata = 9468000; % read all data

time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));

ncfname = [basedir 'TAL_1995_2012_Absolute.nc'];

[x_TAL, y_TAL, z_TAL , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, start_index, ndata, 0);

ncfname = [basedir 'GHC_1995_2012_Absolute.nc'];
[x_GHC, y_GHC, z_GHC , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, start_index, ndata, 0);

ncfname = [basedir 'PBC_1995_2012_Absolute.nc'];
[x_PBC, y_PBC, z_PBC , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, start_index, ndata, 0);

%% find the common data


% L = isnan(x_TAL) | isnan(y_TAL) | isnan(z_TAL) | ...
%     isnan(x_GHC) | isnan(y_GHC) | isnan(z_GHC) | ...
%     isnan(x_PBC) | isnan(y_PBC) | isnan(z_PBC);
% 
% L = ~L;
% 
%% find the common data with GHC and TAL


L = isnan(x_TAL) | isnan(y_TAL) | isnan(z_TAL) | ...
    isnan(x_GHC) | isnan(y_GHC) | isnan(z_GHC);

L = ~L;
%% prepare input data with long continous segments
% idea from http://stackoverflow.com/questions/2212201/


D = diff(L);

b.beg = 1 + find(D == 1);

if L(1)
  b.beg = [1;b.beg];
end

b.end = find(D == -1);

if L(end)
  b.end(end+1) = numel(L);
end

seg_length = b.end - b.beg;

seg_length_sorted = sort(seg_length);
%% get the longest continuous seg




k = find(seg_length == 18124);

% create training data

xTAL = x_TAL(b.beg(k):b.end(k));
yTAL = y_TAL(b.beg(k):b.end(k));
zTAL = z_TAL(b.beg(k):b.end(k));

xGHC = x_GHC(b.beg(k):b.end(k));
yGHC = y_GHC(b.beg(k):b.end(k));
zGHC = z_GHC(b.beg(k):b.end(k));

xPBC = x_PBC(b.beg(k):b.end(k));
yPBC = y_PBC(b.beg(k):b.end(k));
zPBC = z_PBC(b.beg(k):b.end(k));

time_axis = time_array_min(b.beg(k):b.end(k));


% Now make cell structure for ANN input
% Let us assign TAL as local station and GHC as remote station Find local dF, dD

dTAL = 180/pi * atan( yTAL./ xTAL );
fTAL = sqrt(xTAL.^2 + yTAL.^2 + zTAL.^2);

% get remote F
fGHC = sqrt(xGHC.^2 + yGHC.^2 + zGHC.^2);

%get remote delclination
dGHC = 180/pi * atan( yGHC./ xGHC );

% get remote F
fPBC = sqrt(xPBC.^2 + yPBC.^2 + zPBC.^2);

%get remote delclination
dPBC = 180/pi * atan( yPBC./ xPBC );


% make input and output targets. Trend is removed.
% 

% output = detrend(dTAL);
% input = detrend([dGHC zTAL fTAL fGHC zGHC]);

% output = detrend(dTAL);
% input = detrend([dGHC fTAL fGHC]);

output = detrend(dTAL);
input = detrend([dGHC zGHC fGHC]);


% output = detrend(dTAL);
% input = detrend([dGHC fGHC]);


% output = detrend(dTAL);
% input = detrend([dGHC fGHC zGHC dPBC fPBC fTAL]);


%% make validation data set 


k = find(seg_length == 10708);


% create training data

xTAL = x_TAL(b.beg(k):b.end(k));
yTAL = y_TAL(b.beg(k):b.end(k));
zTAL = z_TAL(b.beg(k):b.end(k));

xGHC = x_GHC(b.beg(k):b.end(k));
yGHC = y_GHC(b.beg(k):b.end(k));
zGHC = z_GHC(b.beg(k):b.end(k));

xPBC = x_PBC(b.beg(k):b.end(k));
yPBC = y_PBC(b.beg(k):b.end(k));
zPBC = z_PBC(b.beg(k):b.end(k));



% Now make cell structure for ANN input
% Let us assign TAL as local station and GHC as remote station Find local dF, dD

dTAL = 180/pi * atan( yTAL./ xTAL );
fTAL = sqrt(xTAL.^2 + yTAL.^2 + zTAL.^2);

% get remote F
fGHC = sqrt(xGHC.^2 + yGHC.^2 + zGHC.^2);

%get remote delclination
dGHC = 180/pi * atan( yGHC./ xGHC );


% get remote F
fPBC = sqrt(xPBC.^2 + yPBC.^2 + zPBC.^2);

%get remote delclination
dPBC = 180/pi * atan( yPBC./ xPBC );

time_axis1 = time_array_min(b.beg(k):b.end(k));


% make input and output targets. Trend is removed.


% output1 = detrend(dTAL);
% input1 = detrend([dGHC zTAL fTAL fGHC zGHC]);

% output1 = detrend(dTAL);
% input1 = detrend([dGHC fGHC]);

% output1 = detrend(dTAL);
% input1 = detrend([dGHC fTAL fGHC]);

% output1 = detrend(dTAL);
% input1 = detrend([dGHC fTAL zTAL zGHC fGHC]);

output1 = detrend(dTAL);
input1 = detrend([dGHC zGHC fGHC]);

% output1 = detrend(dTAL);
% input1 = detrend([dGHC fGHC zGHC]);

% output1 = detrend(dTAL);
% input1 = detrend([dGHC fGHC zGHC dPBC fPBC fTAL]);


% plot data

subplot(211);

plot(output,'r');
hold;
plot(input(:,1),'b');
title(sprintf('rms = %5.2f', rms(output-input(:,1))));
subplot(212);
plot(output1,'r');
hold;
plot(input1(:,1),'b');
title(sprintf('rms = %5.2f', rms(output1-input1(:,1))));
